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Abstract. Different options can be used in order to 
measure the shear from observations in the context of 
weak lensing. Here we introduce new methods where the 
isotropy assumption for the distribution of the source 
galaxies is implemented directly on the observed quadrupole| 
moments. A quantitative analysis of the error associated 
with the finite number of source galaxies and with their 
ellipticity distribution is provided, applicable even when 
the shear is not weak. Monte Carlo simulations based on a 
realistic sample of source galaxies show that our procedure 
generally leads to errors « 25% smaller than those asso- 
ciated with the standard method of Kaiser and Squires 
(1993). 
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' 1. Introduction 

One of the most interesting applications of gravitational 
lenses arises in the context of weak distortions of far 
■ sources (e.g., see Tyson et al. 1984, Webster 1985). The 
presence of a gravitational lens (e.g., an intervening clus- 
ter of galaxies) breaks the symmetry of the image of an 
isotropic population of extended sources, with the lensed 
objects preferentially oriented in the tangential direction 
with respect to the center of the lens. This effect is usu- 
ally quantified in terms of the shear introduced by the 
lens. In turn, the shear field measured from the observed 
anisotropy and stretching in a field of distant objects can 
be used to infer the two-dimensional (projected) mass dis- 
tribution of the lens (see Kaiser and Squires 1993; Kaiser 
et al. 1995). 

At present these concepts have found application in 
two different classes of lenses, i.e. galaxies and clusters. 
The work on galaxies has been performed via a statisti- 
cal investigation of a suitable ensemble (Griffiths et al. 
1996, Brainerd et al. 1996), due to the limited number of 
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background sources generally associated with an individ- 
ual galaxy. Relatively nearby clusters, instead, are ideal 
candidates for direct weak lens analyses (e.g., Kneib et al. 
1996; for a review, see Kneib and Soucail 1995; see also 
|Gould and Villumsen 1994); these studies are especially 
important in relation to the problem of dark matter on 
the megaparsec scale, with significant cosmological impli- 
cations. 

In this paper we focus on the measurement of shear 
from observations, as a key step in the process of esti- 
mating the mass of a cluster of galaxies from weak grav- 
itational lensing. We discuss the possible merits of alter- 
native options in handling a given set of data under the 
same conditions for the population of background sources. 
Different methods may introduce different forms of bias, 
depending on whether the data are degraded by seeing 
or by other effects. In this article we give special atten- 
tion to one aspect of the problem which is largely inde- 
pendent of the specific characteristics of the observations 
involved. This important factor in determining the size of 
the error in the shear measurement is the ellipticity dis- 
tribution of the source galaxies. In Sect. 2, with a short 
discussion of the "standard" method of Kaiser and Squires 
(1993), we provide an expression for the expected error in 
a shear measurement as a function of the ellipticity dis- 
tribution of the sources. This result is not restricted to a 
weak shear analysis. Then, in Sect. 3, we introduce a dif- 
ferent method, directly based on the observed quadrupole 
moments, which retains information not only on the shape 
of the observed objects, but also on their angular size, 
and we carry out the related error analysis. The method 
is shown to implement correctly the hypothesis of a ran- 
dom orientation of the sources; it is also shown to be fully 
equivalent to the "standard" method if applied to a popu- 
lation of nearly round source objects. The method is then 
generalized and further improved by introducing suitable 
weight functions to optimize the inversion procedure. In 
Sect. 4 we describe the results of a wide set of Monte Carlo 
simulations applied to realistic distributions of ellipticities 
for the source galaxies. These show that the new method 
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proposed in this paper can be significantly more accurate 
than the standard method in determining the shear from 
the data. The simulations also confirm that the applica- 
bility of the method and the conclusions on its accuracy 
do not depend on the strength of the lens in a relatively 
wide parameter range. 

2. The isotropy hypothesis and the "standard" 
method 

2.1. Notation 

Following standard notation with minor modifications (see 
Schneider et al. 1992), let 9 s be the unlensed position of a 
point source and 9 the position of the corresponding ob- 
served image. We will generally refer to a Cartesian rep- 
resentation of the angle vectors, so that 9 = (61,82). The 
deflection angle (3 is defined through the relation 



9 s = 9-/3(9) 



(1) 



For a single thin lens we can define D os , D^s, and D & as 
the distances between the observer (o), the lens or de- 
flector (d), and the source (s). In this case the deflec- 
tion angle can be expressed as a function of the dimen- 
sionlcss mass distribution k(9) = £(0)/£ c , where £(<?) 
is the projected mass distribution of the deflector and 
S c = c 2 -D s/(47rGDds-Dod) is the critical density; it is 
given by 



(3(9) 



1 f k(9')(9-9') 
\\9-9'F 



2nf 



(2) 



J(9) 



Thus for a small extended source in the vicinity of 9 = 9q 
we have 9 s ~ 9 s (9 ) + J(9 )(9 - 9 ), with 

'86^ 
d9 

l- K (0)+ 7l (0) 72 (9) 

l2 (9) i_ K (0)- 7l (0) 

-<•-«<•» O^" .-%)■ < 3 > 

The complex quantity 7 = 71 + i72 is called shear be- 
cause it is responsible for the distortion of the image. The 
reduced shear parameter is defined as g — g\ + \gi = 
7/(1 — k) and is the quantity actually derived from the 
observations; note that for a weak lens g ~ 7. 

From the surface brightness distribution 1(9) of an ex- 
tended object we can define the total luminosity /, the po- 
sition vector of the center of the image d, and the image 
quadrupole Qij as 



= J 1(9) d 2 6, 



S = j I 91(9) d 2 ' 



Qa = \ J(0 l -S t )(6 3 -5 3 )I(9)d 2 e . 



(4) 
(5) 
(6) 



The quadrupole Qij gives information on angular size, 
shape, and orientation of the observed galaxy. If we are not 
interested in the image size, we may refer to the (complex) 
cllipticity x : 



X = Xi+ 1X2 



Qii-Q 2 2 + 2iQi 2 
Qn + Q22 



x(Q)- 



(7) 



The modulus of \ is related to the shape (|%| = for a 
circular object) while its argument is twice the image ori- 
entation angle (angle between the 6\ axis and the galaxy 
major axis). 

Similar definitions can be given for the unlensed 
source, with 1(9) replaced by the source luminosity flux 
I s (9 s ). For a small object a simple relation holds between 
the observed and the source quadrupole moment: 



s = JQJ , 



(8) 



where we have used the symmetry of J. The corresponding 
relation between \ and x s 1S (here an asterisk denotes 
complex conjugation) 



X 



X + 2g + X *9 2 



x s (x, 9) ■ 



(9) 



f + 2Re( x *g) + \g\ 2 
Thus x s is a function of \ an d 3 only. 

2.2. The standard method (X method) 

The local value of g can be derived from the observation of 
a large number N of galaxies in a small patch of the sky. 
For the purpose, it is usually assumed that the population 
of source objects is isotropic, i.e. that the source galax- 
ies have random orientations. The area of the sky under 
consideration must be small so that the Jacobian matrix 
can be considered to be approximately constant, i.e. the 
same for all the galaxies observed in the area. The source 
galaxies are taken to have £> os /I?ds — constant, which is 
a reasonable assumption for far away sources lensed by a 
nearby cluster. 

The standard method (X method) of Kaiser and 
Squires (1993) implements the isotropy hypothesis by ar- 
guing that the quantity 



X* 



^X> s (x (n) ,s) 



(10) 



vanishes approximately. For a given set of data {x [n) }, 
with n = 1, . . . , N, the condition 



X s -0 



(11) 



is used to determine the reduced shear parameter g. In 
terms of the probability distribution for the source ellip- 
ticities p(x s )> the isotropy hypothesis states that p(x s ) = 
p(|x s |)- Thus we have (x s ) = 0, with a scalar covariance 
matrix (xlXj) = c %- [Here the indices refer to the rep- 
resentation of the complex shear as x = Xi + 1X2 •] Given 
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the range of values for we must have c < 1/2. Then 
the central limit theorem ensures that the variable x s of 
Eq. (10), calculated for N galaxies, must converge to a 
Gaussian with mean and covariance matrix (c/N)Sij in 
the limit N » 1. 

2.3. Error analysis for the standard method 

Here we consider the simple case of a sharp distribution 
of source ellipticities, when the population of objects is 
basically made of nearly round source galaxies. We then 
assume that the probability distribution p is characterized 
by c <C 1. For such distributions, the condition \ s = is 
shown to give a correct determination of the reduced shear 
parameter g (see Appendix A). If we perform a measure- 
ment for a lens with true reduced shear equal to go, then 
the expected statistics for g has mean and covariance given 

by 

(.9) = 90 , (12) 
Cov i3 ( g ) = c(l ~^° |2) 8a. (13) 

Here and in the following, by Cov i:) -(A) = — Xi)(Xj — 
Xj)) we denote the covariance matrix for the random 
variable X with mean X = (X). The expected error, 
cr g = |l — |ffo| 2 | sj c/4JV, is the same on g\ and on g 2 . 
The result stated in Eq. (12) has been obtained also by 
Schneider & Seitz 1995 (see also Schramm & Kayser 1995); 
expression (13) is consistent with the results by Miralda- 
Escude (1991) and by Schneider & Seitz (1995), who, how- 
ever, do not bring out the (l — |#o| 2 ) dependence. 

The error depends on the source distribution via the 
parameter c and shows that the method can lead to an ac- 
curate measurement of g provided N is sufficiently large. 
The error analysis in the general case of a broad distri- 
bution (c < 1/2) is much more difficult. However, if we 
naively extrapolate the conclusions obtained for sharp dis- 
tributions, the largest error on gi and on g 2 should be 
|l — |.go| 2 |/V8A^. In the opposite limit, when c = 0, the 
error vanishes. In fact, when c = all galaxies are circular. 
In this situation the value of g can be derived simply by 
observing the shape of a single galaxy, and thus the error 
on g is zero. 

The dependence of the error on g is quite interest- 
ing. For I go I < 1 the error decreases with increasing \g \, 
and it vanishes when \go\ = 1. This happens because for 
|<7o| = 1 the lens is singular, since det J oc (l — |<7o| 2 ) = 0, 
which means that all galaxies are seen as thin segments, 
with |x| = 1. If the source population is made of nearly 
round galaxies, by observing even a single galaxy with very 
high ellipticity, we can immediately infer that \go\ — 1; of 
course the argument of go is given by the orientation of the 
galaxy. The behavior of the error for \g \ > 1 is explained 
in terms of the go 1— ► local invariance (Schneider & 
Seitz 1995, Seitz & Schneider 1995). Thus it should be 



possible to infer the error for a value of go with modulus 
|<7o| > 1 from the error for g' — l/<?o- I n f &c t (see Ap- 
pendix A), the relative errors for g and g' are the same. 

3. Alternative methods based on the observed 
average quadrupole 

3.1. Simple quadrupoles (Q method) 

We now introduce a new method based on the ob- 
served quadrupole moments rather than on ellipticities (Q 
method). In order to implement the isotropy assumption 
we replace condition (11) by 

/ 1 N \ 

4«E0 (n, ^H- (14) 

^ n— 1 ' 

Note that here the average is performed inside the paren- 
theses, i.e. on the observed quadrupoles. Then it is argued 
that the average quadrupole is traced back to a circular 
source. The method is thus similar to the one used by Bon- 
net & Mellier (1995), but with an important difference in 
that we refer to an average on quadrupole moments, which 
gives different weights to galaxies with different angular 
sizes (see additional comments at the end of Sect. 3.2). 

To see why Eq. (14) holds, we define the probability 
distribution for the source quadrupole moments pq (Q s ) ■ 
Since the quadrupole matrix Q\j is symmetric, pq(Q s ) 
is a function of three real quantities Q\ = Q\i, Q 2 = 
Q12 = Qli' an d 21 = Q22- The isotropy hypothesis states 
that the probability distribution pq is independent of the 
orientation of the galaxy, so that (Q s ) = M Id, with Id 
the 2x2 identity matrix. Furthermore, isotropy requires 
the covariance matrix of Q s to be of the form 

/d + e d- e\ 
Cov fei (Q s ) = e . (15) 
\d- e d + e J 

As the covariance matrix is positive definite, we must have 
d, e > 0. Here the two indices k and I run on the three 
independent components of Q s . It is possible to show that 
e is related to the constant c defined earlier in Sect. 2.2, 
as c = e/M 2 , while d is the variance oi Q\ 1 + Q% 2 , the size 
of the galaxy. 

If we have a large number N of galaxies with quadrupole| 
moments { } , with n = 1,...,N, we can write 

J (^EQ (n) ) J = ^£ gs(n) " MId ' (16) 

^ 71=1 ' ra=l 

where the approximate sign is used because N is finite. 
Let us now calculate the ellipticity \ f° r both sides of this 
equation by applying the definition given by Eq. (7) 

x(j(^EQ ( " } )^x(MId) = 0. (17) 
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Equation (17) is the result stated in Eq. (14), because 

x(JQJ) = x s (Q,g)- 

The reduced shear g calculated through Eq. (14) is a 
good estimate of the true reduced shear g . In Appendix A 
we show that in the case of sharp distributions the quan- 
tity g has mean and covariance matrix given by 



(.9) = 5o , 

n f ^ c ( 1 -|goi 2 )% 
Covy(.g) = — % . 



(18) 
(19) 



Note that these expressions are exactly the same as for the 
standard method. This apparently surprising result can 
be traced to the decoupling of size and ellipticity distri- 
butions for quasi-circular sources. This property is shown 
directly by Eq. (15), where Cov(Q s ) is found to depend 
on two quantities, d being related to the size distribution 
and e being related to the ellipticity. 

3.2. Weighted quadrupoles (W method) 

The simple quadrupole method described in Sect. 3.1 is 
based on the isotropy requirement (Q s ) = M Id, which 
leads to Eq. (16) and thus justifies the use of Eq. (14). 
Since the analysis has shown that g is determined more 
accurately when the ellipticities involved are smaller, we 
may argue that a generalization where a penalty is as- 
signed to galaxies with large |x s | should be able to im- 
prove the method. Thus we introduce a weight function 
W(Q S ) > with the following requirements: 

1. W depends only on the source quadrupoles and is in- 
variant under rotation (consistent with the isotropy of 

pq(Q s ))- 

2. If k is a constant multiplying factor, W(kQ s ) = 
f(k)W(Q s ) (consistent with the fact that the anal- 
ysis of the data is able to constrain J only up to a 
multiplying factor (1 — k)). 

3. W(Q S ) penalizes objects with large |x s |- 

Under these conditions we may now consider as the 
starting point the relation (Q S W(Q S )) = M\y Id, with 
M w a positive number dependent on W, which holds be- 
cause of the isotropy assumption. Then, following the ar- 
gument that has led to Eq. (17), we find 



j(^Y,Q in)w ( J Q {n)j )) J ) -x(m w u) 

^ n=l ' / 



= , (20) 
which indicates that we could replace Eq. (14) by 



/ 1 N \ 

x s (jjY,Q {n)w ( J Q {n)j )>9) = °- 

^ n=l ' 



(21) 



In this method based on weighted quadrupoles (W 
method) one has to proceed by iteration, because the 



quadrupoles JQ^ n \j = Q s ^ that appear in the weight 
functions are not given directly by the data but must be 
guessed first in order to start the iteration procedure. Sev- 
eral alternatives have been considered. One possibility is 
to take the result of the method described in Sect. 3.1 as 
the initial seed for the iteration. Another reasonable choice 
is to take the natural approximation of the weak lensing 
limit, i.e. g ~ —(1/2N) J2 n X^ n \ or tnc estimate obtained 
from g/(l + \g\ 2 ) = -{l/2N)J^ n x {n) applicable beyond 
the weak lensing limit (for a sharp ellipticity distribution; 
see Appendix A). 

As far as the choice of the weight functions is con- 
cerned, we have found that a logarithmic dependence 



W(Q S 



lnf 



X 



(22) 



is simple and suitable for the purpose. 

Note that a weight function could be introduced also 
for the standard X method. As indicated at the beginning 
of this subsection, the proper way to impose the weights 
refers to the source quadrupole moment Q s = JQJ and 
not to the observed quantity Q. Then Eq. (10) would be 
replaced by 



1 N 

-Y, w ( J Q (n)j )x s (x (n \g) ■ 



N 



(23) 



From this point of view, the Q method could thus be seen 
as a special case of weighted X method, with weights 



W(Q S ) = Q s n + Q 



22 ! 



(24) 



as can be checked from the relevant definitions. Similarly, 
the W method can always be transformed into a weighted 
X method, and vice-versa. 

The weighted quadrupole method presents some analo- 
gies with the rejection technique used by Bonnet & Mel- 
lier (1995). Galaxies with large ellipticities are weakly dis- 
torted by the lens (with little information on the shear), 
and contribute significantly to the dispersion. Bonnet & 
Mellier (1995) thus argue that only galaxies with small 
ellipticities should be retained and galaxies with elliptic- 
ities larger than a given threshold value Xcut should be 
rejected. Hence rejection can be considered as a particular 
weighted method with a step-like weight function. How- 
ever, the rejection apparently used in the above-mentioned 
article refers to the observed (rather than the source) ellip- 
ticities. This in turn introduces an undesired anisotropy in 
the galaxy sample, and thus a bias in the measured shear. 
This point is easily clarified in the weak lensing approxi- 
mation. In this case Eq. (9) becomes simply x s — X + 270 
so that the lens acts as a translation on the ellipticity. A 
rejection on |x| then becomes a rejection on |x s — 270 1, i.e. 
a rejection of all galaxies that in the x s complex plane are 
outside a circle of radius x cu t centered on 270. Therefore, 
using the X method, the expected value for 7 is reduced, 
i.e. 



(7) 



70 - ^P(Xcut)Xcut70 



(25) 
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where T is the fraction of unrejected galaxies and Xcut 
is the maximum allowed value for |x|. This expression is 
valid for Xcut < 1 — 270 and to first order in 70. Simi- 
lar comments could be made on the rejection using the Q 
method, but here (7) depends in a more complicated way 
on the probability distribution pq. In conclusion a rejec- 
tion referred to the observed ellipticities introduces a bias 
on 7 in the direction of smaller values. 

One might also try to relate the W method to maximum-| 
likelihood methods used in this context (e.g., see Bartel- 
mann et al. 1996). The relative merits may depend on the 
characteristics of the population of source galaxies (see 
also Sect. 4). 

The arguments discussed after Eq. (22) have found 
convincing demonstration in the simulations outlined be- 
low in Sect. 4. 

3.3. Advantages with respect to the standard method 

The following is a short qualitative discussion of the three 
methods introduced above in relation to the properties of 
the background source galaxy population. 

Our starting point is the fact that the expected error 
on g is proportional to \/c, i.e. it is larger for broader dis- 
tributions. This suggests that the error is dominated by 
the contribution of relatively flat source galaxies. This has 
led us to introduce the weighted quadrupole method. On 
the other hand, the simple quadrupole method should also 
perform better than the standard method. [We should em- 
phasize (see comment after Eq. (19)) that the apparent full 
equivalence of the X and Q methods suggested by Eqs. (12, 
13) and (18, 19) holds only in the limit c — > 0, where 
size and ellipticity decouple from each other, as shown by 
Eq. (15).] This may be argued in the following way. Con- 
sider a simple case where all source galaxies are flat disks 
with the same luminosity and the same size. Then, for 
a random spatial orientation, all galaxies will appear as 
ellipses, with identical major axes and minor axes rang- 
ing from to the size of the major axes. Therefore flatter 
source objects occupy smaller areas of the sky. In contrast 
with the standard method, the quadrupole method takes 
the galaxy size into account, which should lead to more 
accurate results. The argument can be extended to the 
case where different luminosities and sizes are involved. 
Furthermore, it should also be applicable to a more re- 
alistic galaxy population, since galaxy fields are usually 
dominated by spiral galaxies. This will be confirmed by 
our simulations (see Sect. 4). 

It is also interesting to consider the case of a source 
population made of spheroids (i.e. axisymmetric ellip- 
soids). For a population of oblate galaxies we expect a 
behavior similar to that of disk galaxies. Thus if elliptical 
galaxies are generically oblate, their contribution would 
not alter the argument in favor of the quadrupole method. 
The situation is completely different for prolate spheroids, 
which behave in the opposite way; for a population domi- 



nated by prolate objects better results would be expected 
from application of the standard method. 

Finally, it would be important to assess the error in- 
volved in the measurement of the ellipticity, but this would 
depend on a number of conditions characterizing the spe- 
cific set of observations under investigation. Since the 
scale length of the lensed galaxies is often smaller than 
1", seeing can be one important factor for ground ob- 
servations. Seeing makes galaxies rounder and thus leads 
to an underestimate of the shear. Even if seeing effects 
can be partially resolved by means of special algorithms, 
such as maximum-likelihood and maximum-entropy im- 
age restoration (Lucy 1994; see also Bonnet & Mellier 
1995, Kaiser et al. 1995, Luppino & Kaiser 1997, Villum- 
sen 1995), the error on the measured ellipticity should be 
larger for smaller galaxies. This is again in favor of the 
use of the quadrupole method, which downplays the role 
of small galaxies. For simplicity, seeing and other sources 
of error, such as Poisson noise, sky luminosity, and pixel- 
ing are not considered in this paper and in the simulations 
described below. 



3.4- Weight optimization 

The logarithmic choice of Eq. (22) is one simple option 
compatible with the requirements listed at the beginning 
of Sect. 3.2. One may wonder whether it is possible to con- 
struct an optimal weight for a given population of source 
galaxies. Suppose that the probability distribution pq (Q s ) 
is known. Then, at least in principle, one could try to find 
the weight function that minimizes the error on g. 

This interesting possibility is clarified by the follow- 
ing example. Consider the weighted X method described 
by Eq. (23) with a weight function depending only on 
the determinant D s = detQ s . In this case the second re- 
quirement of Sect. 3.2, i.e. W(kQ s ) = f(k)W(Q s ) just 
demands that VF(D S ) be a homogeneous function. There- 
fore W can be taken of the form W(D S ) oc [D s ) n . Let us 
now introduce the distribution probability pd(x s ,D s ) = 
Pd(\x s \j D s ) for the source galaxies. Obviously this distri- 
bution is related (in a complicated way) to the quadrupole 
distribution pq. Using the relation D s — (det J) 2 D, where 
D = det Q, and the form of W(D S ) , the variance of g can 
be shown to be proportional to the variance of x s divided 
by the square of W^D S ) , where D s is the mean value of 
D s (see Eq. (23) and the comment after Eq. (A7), B being 
proportional to W(D*)). In other words the quantity to 
be minimized is 

Covg) _ x 
[W(D^)] 2 N[W(D^)] 2 

x J[W(D s )] 2 \ X s \ 3 p D (\ X s \,D s )d\x s \dD s . (26) 
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Recalling that W(D S ) cx {D s ) v , we can find the value of 77 
that minimizes Eq. (26) provided that we know the source 
distribution p D . 

The value of r\ that minimizes expression (26) can be 
easily obtained numerically. An approximate solution can 
be found by expanding D s near its mean value D s . To 
second order 



1 D S {(D 

Wheat = 7 ~ 



W)\ x s \ 2 ) 



2((D S - D s ) \x s \ 2 ) 



(27) 



This should lead to the most accurate determination of g. 
The second term of Eq. (27), with its minus sign, takes 
into account the correlation between the galaxy size D s , 
and its ellipticity x s - If, as in our simulations, the largest 
galaxies are the ones with smallest ellipticity, then this 
term is positive. In fact, here it is convenient to use a large 
value for r\ in order to penalize smaller galaxies, which 
should be the ones with a large ellipticity. Curiously, the 
case where size and ellipticity are fully uncorrelated in the 
source population would still require a non-trivial weight 
function (r/bost = 1/4). 

Of course the weight optimization should be performed 
in both the source size (D s ) and ellipticity (x s ). This task 
is, obviously, much more difficult, mainly because of the 
non-trivial dependence of x on % s , given by Eq. (A8). 

4. Monte Carlo simulations 

4-1- Source galaxies 

To test the various methods we have performed a series of 
Monte Carlo simulations of measurements of the reduced 
shear parameter g, by generating N galaxies as a realiza- 
tion of two types of source galaxies. At this stage it is not 
completely clear what is to be considered a realistic popu- 
lation of source galaxies. Below we focus on natural cases 
that we could draw from the literature. 

The first galaxy distribution (Population A) is char- 
acterized by a source ellipticity distribution inferred from 
about 6000 galaxies observed in a single frame obtained 
at the 200-inch Hale telescope at Palomar (Brainerd et 
al. 1996). Calling q the axis ratio of the galaxy (with 
< q < 1), the empirical probability distribution for q 
is 



p q (q) ~ 64 g exp(-8<?) . 



(28) 



The data by Oderwahn et al. (1997) basically confirm the 
plausibility of this choice. The associated distribution for 
the complex ellipticity x s is shown in Fig. 1. The source 
galaxies are largely dominated by disk galaxies, with sig- 
nificant contributions from ellipticals. In order to spec- 
ify the size distribution, for simplicity we have followed 
Wilson et al. (1996) in adopting exponential luminosity 
profiles for all the source objects, with scale-length h uni- 
formly distributed in the range [0.25, 0.65] arc-seconds. Fi- 
nally, the position angle has been chosen randomly with 



Fig. 1. The probability distribution of the source ellipticity 
for the galaxy populations A and B described in the text. The 
distribution is isotropic, i.e., p(x s ) ~ p(|x s |)- Note that in 
our notation the distribution is normalized with the condition 
2 7 r/p(|x s |)|x s |d| X 8 | = l. 

uniform distribution. Wc have assumed that the three dis- 
tributions (ellipticity, size, and position angle) are inde- 
pendent. A correlation between ellipticity and size might 
be present in a more realistic population of sources. Based 
on this choice of p q (q) we have c ~ 0.0606. Thus the 
expected error in the determination of g for both the 
standard method and the simple quadrupole method is 
« 0.015/AT. 

The second simple population of source galaxies (Pop- 
ulation B) has been generated by assuming that all galax- 
ies are flat disks of the same size and with random orienta- 
tion (see also Bonnet & Mellier 1995). A straightforward 
calculation shows that the probability distribution for q is 
uniform, i.e. 



Pq(Q) = 1 



(29) 



In this case we have c ~ 0.2146. The associated distri- 
bution for x s (see Fig. 1) shows a pronounced peak near 
|x| = 1. This distribution is not sharp in the sense of 
Sect. 2.3, because of its relatively large value of c. 

4-2. Simulations and results 

The lens properties arc then specified by means of the 
value of go] the value of (1 — k) does not affect the results. 
For a given lens, the observed quadrupole moment 
of each source is calculated by inverting Eq. (8) 



(30) 



From the set of observed quadrupole moments { } 
the methods X, Q, and W described in Sect. 2 and Sect. 3 
allow us to determine the measured shear g. All methods 
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Fig. 2. The graphs show, for population A, the mean error on 
g as a function of go- The error has been obtained from 10, 000 
simulations (N = 16 galaxies). Frame (a) shows the results 
for the standard method. The three methods are compared in 
frame (b). The results agree very well with the |l — \go\ 2 \ law. 

involve the resolution of an implicit equation (Eq. (11), 
Eq. (14), and Eq. (21)): this has been done with a simple 
Newton-Raphson algorithm (see Press et al. 1992). For the 
results shown below the W method has been implemented 
in the form of Eq. (22); tests on other types of weight 
functions, not shown here, have also been performed. 

The entire process has been repeated a large number 
of times (typically 10,000), each time based on a differ- 
ent realization of the N source galaxies. Thus for each 
method, the mean of the measures of g and the related 
errors (covariance matrix) have been calculated. The sim- 
ulations immediately provide a simple check on the mean 
of x s and Q s and on their covariance matrices, consistent 
with the relation c = e/M 2 stated in Sect. 3.1. Then we 
have compared the average errors of each method with the 
expected error given by Eqs. (13) and (19). 

The main results for population A of galaxies are 
shown in Fig. 2, confirming the anticipated dependence 
on |<7o|; the expected diagonal character of Gov ij(g) for all 
the methods is also recovered. 

As we can see from Fig. 3, the errors for the standard 
method (X) follow the 1/y/N law. They agree quite well 
with Eqs. (13) and (19) for all go values. From these re- 



Fig. 3. The mean error on g versus the number N of galax- 
ies used in the reconstruction process. Symbols are the results 
of the average error obtained with 10, 000 simulations for the 
three specified methods, while the line is the l/y/N law given 
by Eqs. (13) and (19). The value of g used is g = 0.2 + 0.2i. 
The trends observed at N ~ 30 have been checked to persist 
well beyond the range of N shown here. 

suits on the standard method, the galaxy distribution A 
is found to be sufficiently sharp. The error for the sim- 
ple quadrupole method (Q) is instead about 15% lower 
than in the standard method for every value of g and N 
(see the short discussion given in Sect. 3.3). The weighted 
quadrupole method (W) shows even smaller errors: these 
are about 28% smaller than that of the standard method. 
The fact that the X symbols are closest to the theoretical 
line just demonstrates that the asymptotic approach at 
the basis of Eqs. (13) and (19) is more "natural" for the 
X method. 

The results based on population B are qualitatively 
similar. The relevant ellipticity distribution is rather 
broad, and thus we expect Eq. (13) to give only approx- 
imate results. In fact, actual errors are about 30% larger 
than those stated in Eq. (13). For this galaxy population 
the differences among the three methods are more signif- 
icant. The error in the simple quadrupole method is 27% 
smaller than that in the standard method, while in the 
weighted quadrupole method it is 38% smaller. 

In order to further test the sensitivity of our results to 
the properties of the source galaxy distribution, we have 
run simulations on another class of populations, not nec- 
essarily following empirical suggestions (as was the case of 
Eq. (28)), but rather mathematical simplicity. This class 
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of populations has been produced by imposing a Gaussian 
distribution for x s (suitably truncated at |x s | = 1). One 
case, that we have called Population C, is characterized 
by a truncated Gaussian with covariance c = 0.0606 and 
by the same size distribution as for population A. Here 
we find that the Q method performs better than the X 
method, with error smaller by w 20%; curiously, the W 
method implemented with Eq. (22) here leads to errors 
comparable to those of the X method. It thus appears 
that the optimal choice of reconstruction method may be 
sensitive to the detailed behavior of p q (q) around q = 1. 

The specific examples addressed by the simulations can 
also be interpreted in terms of an "effective population" of 
source galaxies, determined by the adopted weights (which 
are a function of x s )- In other words, good weights act as 
if the dispersion in the distribution of x s were reduced. 
Of course this is only an intuitive way of describing how 
the weights operate in the process, since the resulting "ef- 
fective population" depends on the unknown shear. In ad- 
dition, one should also introduce an "effective N" to de- 
scribe the impact of the weights not directly related to the 
shear. 



5. Conclusions 

In this paper we have analyzed in detail different meth- 
ods aimed at obtaining the reduced shear g from the ob- 
served ellipticities of a set of galaxies and we have focused 
on the issue of the accuracy of such measurement. Monte 
Carlo simulations have demonstrated that realistic distri- 
butions in general favor a new method introduced in this 
paper, called the quadrupole method. In particular, the 
weighted quadrupole method has been shown to perform 
best (the gain is of the order of ss 15-30% for the simple 
quadrupole method and can be of the order of w 30-40% 
for the weighted quadrupole method). This in turn leads 
to an improvement on the mass determinations, which is 
anticipated to be by a similar factor. 
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helpful discussions and suggestions. This work has been 
partially supported by MURST and by ASI of Italy. 



Appendix A: The expected covariance matrix for 
a sharp ellipticity distribution 

In order to calculate the covariance expected in a mea- 
surement of g with the methods described in the main 
text, we use a general expression known in the theory of 
errors (see Taylor 1982). Let A be a multidimensional 
random variable and let Y = ,f{X). Suppose that / is 
a smooth function and that the probability distribution 
of X is sharp, so that its covariance matrix Cov(A) is 



small. Then we have (basically this is an application of 
the saddle-point integration method) 



(Y)cf((X)), 
Cov(F) ~ CCov(A)C T 



(Al) 
(A2) 



Here C is the Jacobian matrix for / calculated for the 
mean value of X, i.e. 



C 



df_ 

dX 



(A3) 



X = {X) 



If Y is implicitly defined by F(X, Y) = 0, we can equally 
well apply the previous relations using the inverse function 
theorem. In this case the mean value of Y is implicitly 
provided by F((X), (Y)) = 0, while C = -B^A with 



OF 

A= ox 

B= 9 -l 
dY 



X = (X),Y={Y) 



X = (X),Y=(Y) 



In conclusion in this case we have 
F((X),(Y}) ~0, 

Cov(F) ~ (B- 1 A)Cov(X)(B- 1 A) T 



(A4) 
(A5) 

(A6) 
(A7) 



If we take Y = (Y) = constant and A to be a random 
variable, then we have Cov(F) = ACov(X)A T . There- 
fore, from Eq. (A7) we may recognize that Cov(F) ~ 
B- 1 Cov(F)(B- 1 ) T . 

Based on these expressions we can easily prove Eqs. (12,| 
13) for the standard method, and (18, 19) for the sim- 
ple quadrupole method. For example, for the standard 
method we will identify Y with the reduced shear and 
F with x s 7 under the condition that the lens be charac- 
terized by true reduced shear, Y = g - 

A.l. Standard method 

In the standard method the reduced shear g is calculated 
from Eq. (11). In principle this equation depends on N 
complex random variables x^ an d on the unknown re- 
duced shear g. The expected average (%) for the observed 
ellipticity of every galaxy x^ is easily obtained from the 
inverse of Eq. (9), that is 



X = 



X s - 2g + x s * g 2 



1 - 2Re(x s *.g) + \g\ 2 ' 
By setting x s = (x s ) = 0, this gives the mean 

<%>=- ^ 



1 



l.9o | 



(A8) 



(A9) 



As in the main text, g represents the real value of g. We 
stress again that this equation is valid only in the case of a 
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sharp source ellipticity distribution (c«l; see Schneider 
& Seitz 1995 for the relevant analysis). 

The expected mean of measurements of g obtained 
from Eq. (11) obeys the relation 



1 N 

^£x s «x>,<s»=o. 



(A10) 



n=l 



Using the calculated mean (x) of Eq. (A9) above we thus 
prove Eq. (12). 

Following Eq. (A7), the covariance matrix for g can be 
written as 



N 



Cov(g) = J2( B ~ lA ") Coy W( B ~ lA ") T ■ 



(All) 



n=l 



The matrices A n and B are the partial derivatives of 
X s (X, 9)- 



An. — 



dx s ({x (m) },g) 



i d x s (x {n \g) 



X (m) =(x),S=(9> 



B 



N dx (n) 

dx B ({x {m) },g) 



dg 

dx s (x,g) 



x (n) = ( x ) lff =(g> 



x=(x),a=(a) 



(A12) 



Dg 



X=(x},S=(ff) 



(A13) 



Note that the matrices A n = (1/N)A are all equal. The 
covariance matrix Cov(x) can be expressed, by application 
of the inverse function theorem, as 



Cov( x ) = A- 1 Cov( x s )(A- 1 ) T = cA-^A-'f 
so that Eq. (All) becomes 

Cov( 5 ) = ^S- 1 (B- 1 ) T . 



A simple algebraic calculation leads to 
2 



B 



Id 



(A14) 



(A15) 



(A16) 



1 - M 2 
and thus we recover Eq. (13). 

A. 2. Simple quadrupole method 

To show that the simple quadrupole method leads to sim- 
ilar results, we start by considering the expected distribu- 
tion for the observed quadrupoles Q. From Eq. (16) in the 
sharp distribution limit we have 



for every galaxy. Here J" 2 indicates (J 1 ) , i.e. the square 
of the inverse matrix of J. The expected average of g is 
calculated by solving Eq. (14) with g replaced by (g), i.e. 



X s ((0), (.9)) = 0. 



(A18) 



Since x s (Q,.g) = x{JQJ)i with J any matrix with associ- 
ated reduced shear equal to g, we recover Eq. (18). 

In order to calculate the covariance matrix of g, we use 
Eq. (A7). From Eq. (A2), the covariance of Q (recall the 
definition Qij = Qi+j-i) can be written as 



Cov(Q) = CCov(Q s )C T , 
where C is the Jacobian matrix 
dQ(Q s ,g) 



C = 



(A19) 



(A20) 



We observe that the argument of x s used in Eq. (14) is 
Q = (1/N) J2 n Q {n): this allows us to usc Q as X variable 
in the analysis. Thus the matrices A and B of Eqs. (A4) 
and (A5) are given by 



A = 



B 



dx s (Q,g) 



dQ 

dx s (Q,g) 



Dg 



Q=(S),s=(s> 



Q=(Q),g=(g) 



(A21) 
(A22) 



where the average values of Q and g are stated in 
Eqs. (A17) and (18). Note that the matrix B defined here 
is the same as the one used for the standard method and 
defined in Eq. (A13). The covariance matrix for g is then 



Cov( 3 ) = (B-'A)Cov(j- £ Q(»>) (B-'A) 

^ 71=1 ' 

= ±(B- 1 A)CCov(Q s )C T (B- 1 A) T 



Therefore we have 

dx s dQ dx s 



AC = 



dQ dQ s 
and thus 

ACCov(Q s )(AC) T = Cov(x s ) = c Id . 
Using Eq. (A16) we finally obtain Eq. (19). 



(A23) 



(A24) 



(A25) 



(Q) = MJ- 



(A17) 
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A. 3. The dependence on g 

From local observations it is not possible to choose be- 
tween two different solutions for the reduced shear related 
by g' — l/g*. We thus expect that the expression for the 
error will incorporate this invariancc property. 

The g i ^ l/g* transformation is equivalent to a change 
of sign in one eigenvalue of J (sec Schneider & Seitz 1995), 
which is not observable. Every reconstruction method is 
then bound to give two solutions for the reduced shear. 
As g' is related to <?, we can calculate the expected error 
on g' given the error on g: 



Cov( ff ') 



dg 



Cov(g) 



dg 



A simple calculation gives 



dl\ = J_(g!-gi -29192 
dg) M 4 ^ -25132 9l-9i 

so that, as Cov(g) oc Id, we have 



(A26) 



Cov(ff') - T-TyCovfo) 



c (i/i. 9 i 2 -iy 



AN 



Id 



c 1 - Iff 



l\2\ 



AN 



Id . 



(A27) 



Hence our expression for the error on g reflects the g i— » 
l/g* invariance. 

Appendix B: Basic structure of the simulation al- 
gorithm 

The simulation algorithm, briefly described in Sect. 4, is 
composed of three different parts. 

For a given galaxy distribution, we have generated N 
source galaxies (quadrupoles) using the transformation 
method (see Press et al. 1992); the rejection method is less 
convenient here. Then we have calculated, using Eq. (30), 
the observed quadrupoles and cllipticities. 

Different methods (X, Q, and W) have been applied by 
solving respectively Eqs. (11), (14), and (21) with a sim- 
ple Ncwton-Raphson algorithm. This algorithm requires 
an initial point for the unknown g. This has been chosen 
using Eq. (A9). Through the entire simulation algorithm 
we have supposed to be able to distinguish between the 
two solutions g and l/g* . In practical cases this is easy if 
the lens is non-critical. For critical lenses this ambiguity 
can be resolved only globally. Note that this assumption 
has been implicitly made throughout the paper (e.g., con- 
sider the justification of Eqs. (A10) and (A18)). 



The previous steps have been repeated a large number 
of times with the same go. For each individual "simula- 
tion" we have memorized the calculated value of g. Finally, 
the mean and the covariance matrix have been calculated 
from the data base of all the available results. 

The algorithm has been implemented as a C++ code 
running on a IBM RISC System/6000 590 machine. 
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